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Abstract 

We derived a number of numerical methods to treat biomolecular systems with multiple time 
scales. Based on the splitting of the operators associated with the slow-varying and fast-varying 
forces, new multiple time-stepping (MTS) methods are obtained by eliminating the dominant terms 
in the error. These new methods can be viewed as a generalization of the impulse method 1, 2]. 
In the implementation of these methods, the long-range forces only need to be computed on the 
slow time scale, which reduces the computational cost considerably. Preliminary analysis for the 
energy conservation property is provided. 
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I. INTRODUCTION 


Simulations based on dynamics models for bio-molecular systems have become standard 
computational tools for studying small-scale motions or calibrating coarse-grained models. 
A core element in these simulations is the numerical integrator for solving the underlying 
differential equations. A well known challenge is the presence of multiple time scales, usually 
arising from the various types of molecular interactions. For example, the bond length and 
bond angle contributions constitute the most dominant terms, and they determine the fastest 
time scales. As a result, the step size of a numerical integrator has to be selected accordingly. 
For instance, the step size for the Verlet’s method needs to satisfy (sj], 

St < —, (1) 

^max 


where a; max is the maximum frequency. The bound on the step size is typically very small 
(femto-seconds or less). 

In practice, the small step size St imposes a significant limitation on how long the sim¬ 
ulations can be conducted. In particular, the most expensive part of the computation is 
the force calculation. While interactions via changes of bond length and bond angles are 
short-ranged and easy to compute, there are long-range interactions (electrostatic) that take 
up considerable CPU time. 

A remarkable approach to overcome this difficulty is the multiple time stepping (MTS) 


method, referred to as Verlct-I in [1] and r-RESPA in [2[] . In sharp contrast to conventional 
integrators, this method involves multiple step sizes, which correspond to the time scales that 
forces of different nature determine. The MTS method is also known as the impulse method. 
The implementation of the impulse method is quite simple: For each small time step St, one 
integrates the ODEs with only the fast force, and for each large time step At, one updates 
the momentum by applying the slow force, which is often described as a half-step ‘kicking’, 
‘oscillating’, and another half-step ‘kicking’. The method can be mathematically written as 


Q. 


a symmetric splitting of the Liouvi 
software packages (e.g., TINKER 


le operator jsj], and it has been implemented in several 
4]). From a practical viewpoint, the impulse method 
eliminates the need to evaluate the long-range forces at every step with step size St, which is 
the most expensive part of the dynamics simulations. As a result, the computation is sped up 
considerably. Another important mathematical property is the symplectic structure, which 
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for Hamiltonian systems, is critical to preserve the total energy and produce the correct 
statistics. 


ica 


issues. For in- 


1 , 6]. Furthermore, 


Despite the popularity of the impulse method, these are known pract 
stance, instability has been observed for some particular choices of At jtj. 1 
the order of accuracy is largely unknown. This has motivated several generalization and im¬ 
provement of the impulse method. For instance, it is generalized to incorporate molecular 
interactions of multiple scales (> 3) by Procacci and Marchi [3], in which the evolution 
operator is split into operators that represent the forces of different magnitude. The LN 


method 


8] combines the slow and fast forces via extrapolation, and Langevin dynamics that 


represents a heat bath is introduced to reduce energy drift. Another important development 


is the MOLLY mot hod 


3 , 


TO], where the impulse method is modified by properly averaging 


the slow force. The goal of the LN and MOLLY methods has been to overcome the resonance 
instability of the impulse method, which arises when the slow time step At is a multiple of 
half of the period of the fast mode. 


This work continues with the operator-splitting procedure that was used to derive the 
original impulse method 3- More specifically, we seek splitting methods that involve more 
fractional steps. In order to determine the coefficients of the splitting methods, a high-order 
error expansion is needed, which in general, is a tremendous challenge, especially when 
multiple fractional steps are involved. Based on the symbolic code that we had recently 
developed jll|, we are able to obtain the error up to any prescribed order. This constitutes 
the basis for determining the coefficients to maximize the order of accuracy. 


The different magnitude of the fast and slow forces has motivated us to introduce another 
parameter e, which indicates the resulting time scale separation. Our analysis reveals that 
some of the terms in the error are weighted by large factors such as 4j. This prompts us to 
re-consider their role in the error. By eliminating dominating terms in the error, we obtain 
three new impulse methods, which involve two to four fractional steps. Numerical tests 
suggest that the new methods have improved accuracy. Furthermore, these methods also 
have better performance in conserving the total energy. These numerical observations can 
be interpreted with some preliminary analysis. 
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II. MATHEMATICAL FORMULATIONS 


In this section, we present the mathematical framework for deriving new impulse methods. 


A. Problem setup 


To explicitly indicate the fast-varying force in the dynamics, we consider a dynamical 
system in the following form, 


Xi = Vi , 


Vi = + fi{x), i = 1,2, • • • , N. 


( 2 ) 


In particular, the fast time scale is indicated by the parameter e, with e « 1. 

Numerical methods for such ODEs, especially the one-step methods, are based on the 
approximation of the evolution operator, which for this autonomous system can be expressed 
as, 

(x(t),v(t)) =e tSf (x,v). (3) 

Here (x, v ) denotes the solution at t — 0. Furthermore, the differential operator is defined 
as, 


& = 


( 4 ) 


All the partial derivatives are defined with respect to the initial data (x,v). 

An important class of methods have been based on the splitting of the operator Jzf in the 
semi-group operator (|3|h As a one-step method, it suffices to consider the approximation 
formulas for e AtJf , since the same formula can be applied to all the following steps. For the 
present problem, we define 

1 

C/J C/)X . C/)v 


=Y. v ‘ d « 

i 

i 

3 

The popular impulse method is based on the following specific splitting scheme 


( 5 ) 


Q, 


= AtS? ^ TAiif 2o A^ijA(y 2 


ea 


l e2 


( 6 ) 
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In this paper, we seek splitting methods of the general form, 


k 

Atse _ e C’ At ^ e (L ' Am 

i =1 


Here the right hand side is interpreted as, 


( 7 ) 


e d At^ e di ALA e c 2 e d 2 Am _ _ _ e c k AtJ? 2 e d k Am 


For example, in the impulse method ([6]), we have C\ — c 2 — |, d\ — 1 and d 2 = 0. 

Notice that the step e CiAt can be implemented exactly, since the corresponding differ¬ 
ential equations have explicit solutions. On the other hand, the step e c,A< -^ 1 requires further 
approximation. The corresponding differential equations are, 


Xi = Vi , 

Vi = \gi(x) 


( 8 ) 


For instance, it can be approximated by the Verlet’s method with step size St -C e. Namely, 


CiAtA'i 


St c/>v £+ c/>x dt c/> 

A C.C' 1 (J t'oC' "1 _ A oi' 1 


St , 


(9) 


in which ht = In this case, the error is expected to be on the order of 0(St 2 ), which is 
negligible compared to the larger time step At. Therefore, we will simply assume that the 
step e c<At5fl is implemented exactly, and focus primarily on the error from the splitting at 
the larger steps At. 


As discussed in the introduction, the advantage of introducing such splitting methods 
arises when the fast force g is short-ranged, the computation of which is much less expensive 
than that of f(x). Therefore, compared to a direct discretization of the ODEs ([2]), where 
f(x) is computed at every step, the splitting methods can be implemented at a reduced cost. 
Furthermore, one can show that all such splitting methods are symplectic, when the ODEs 

he computer implementation of these methods is 


form a Hamiltonian system 


12], Finally, t 


quite straightforward, as demonstrated in 


13 


14]- 


On the other hand, the accuracy of this method has not been extensively studied. Our 
mathematical formulation relies on an explicit representation of the error for the splitting 
method (J7|). More specifically, given the two differential operators and the number 
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of fractional steps k, and the coefficients (ci, C 2 , • • • Ck) and (di, d 2 , ■ ■ ■ , dk), we seek SB, such 
that: 

e c . At ^ e t:Am = • • • Cfc ; d 0 , ■ ■ ■ d k )e At ^ (10) 

1=1 

where the near-identity operator SB can be expanded as follows, 

^(co, • • • c-k] do,--- dk) ~ y = X x A t + X 2 X t 2 + ■■■ XAt n + ■■■ ■ (11) 

Here A is the identity operator, and Xj contains terms with j-folded operator multiplica¬ 
tions, with coefficients represented by [c\, ■ ■ ■ Ck, d\, ■ ■ ■ ,dk), he., cxOid^X^XiX^Xi £ X4. 
All terms appeared in Xj are linearly independent. With the same effort, one can also 
express (ITTh in the following form, 

JJ e CtAt ' % e dlAm = e AtS/ 'M{co, ■ ■ ■ c k ; d 0 , ■ ■ ■ d k ). (12) 

1=1 

The operator SB can be expanded in a similar format. 

The local accuracy of the splitting method is determined by the magnitude of the first few 
none-zero operator coefficients on the right hand side of (TlTh . Therefore, it is necessary to 
derive the explicit expressions for the first few coefficients. This procedure will be described 
in the next section. 


B. Finding explicit forms of SB 

Without loss of generality, let us consider two operators srf and SB. The starting point 
of our analysis is the following expansion, 

srf _ (j7> ( 1 o\ 

6 6 •> \ 

or alternatively, 

6?) _ „S8 0 —s£—£8 (-\A\ 

— e e e . (14) 

We want to express this approximation in operator multiplication form, instead of the 
exponential form in [lfi]], because multiplication form often gives back cleaner expressions 
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ll|. Direct computation shows that: 



{siSB - SBsrf) 

- Cs/ 2 .% - Xy/tA.i/ - e/SS 2 + 2SSe/ 2 + 2SgsASS - lA-.i/) 


( 15 ) 


More terms are available, simply by Taylor expansion of each term. 

We now consider further splitting of sS + SB, in the general form of, 

k 

]Je c ^ e di ®. ( 16 ) 

i =1 

The key to obtain the expansion of the error for (11611 is to repeatedly use the above 
formula: 

k 

e e 

i =1 


=&{curf,<h3g)&{cutf+di3g,C2stf) eClS, * +dl3g+C2 ^ I! eC ^ edi * 


_ C/7) crt> „cisrf+C2£f+di3g+d,23B I I n diS§ 

,d\38)’y4'(c\srf+d\3B,C2srf) t ^'{c\£/+C2S?+d\3B,d238)& XX ^ ^ 

i =3 


k 

= J^[ &(C i ii/+D i - 1 &,di£g)&(Ci£/+Di&,ci +1 &)e Ck * /+Dk ^ 

2=1 

^( Cl ,---c k -d u --- ,d k )e c ^ +D ^ (17) 

where 

i i 

Ci = y ^ Cj , d, = y ^ dj , 

j =1 j =1 

for 1 < i < k, are the partial sums of the coefficients, and we define Co = D 0 = 0. 

This systematic procedure will be applied to analyze the splitting methods (El) for the 
multiscale ODEs (J2]). In particular, we let sf = AtJ*f 2 and SB = AtJSi. We require that 
Ck = Dk = 1, which is the consistency condition for one-step methods Q . This leads to 
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the error for (fTOl) and (fill) . In particular, we now have, 


^(ci, • ■ -c fc ;rfi, • • • ,4) 



fc 


( 18 ) 


Z=1 


Clearly, it would be a lengthy procedure to carry out the multiplication of these operators. 


Fortunately, we have developed a symbolic code jll] to obtain the expansion of the error. 
This constitutes the basis to examine the coefficients of the error, from which the order of 
the accuracy can be determined, controlled and improved. In the next section, we discuss 
numerous cases. 

C. The selection of the coefficients a and di based on the error expansion 

We now return to integrators with the general form (ITjh We will examine the cases k = 2, 
k = 3 and k = 4 separately. 

1. Error expansion for the case k = 2. 

We first consider the case k = 2. Based on the analysis from the previous section, we 
found, 


&{c h C 2 = 1 — ci, 4,4 = l - 4) - 0 


((Cl - 1)4 + - AJA 2 )At 2 


(19) 


+ ( ([ ' { \ Cl) - y)(J^J4 - t 3 

2 o 


+ ( d ' (C ' 2 1] + i)(2^24A? 2 - 2Sf 2 ^ 2 + Af 2 Af 2 Afi)Af 3 + ... 


In the traditional impulse method, C\ — c 2 = and it may appear as if the error is 
0(At 3 ). But equation (J5J) suggests that is which clearly indicates that a further 

inspection of the order of the error terms is needed. To this end, let us define, 


0 21 A24J4-J4J4, 

03! AA\A\/A 2 - Z.A\.A 2 .A\ + A 2 A\A, 
03 2 =A X A 2 A 2 - + JA 2 fA 2 fA i. 


( 20 ) 


We first begin with the observation that, 




Lemma II. 1. The following identity holds, 


= S^Sf 2 . (21) 

The equality can be checked directly, 

2 ^2 ^ - 25^2% = EDaa-w.) 

* 3 

By assuming that the flow of the ODEs is sufficiently smooth so that df. v . 
get S&Sf? = 2 ^ 2 %. 

Following this calculation, we find that, 

Theorem II.2. The order of the operators in (12U1) is given as follows, 

1)®21 = 0(1), 

2) 3*31 = 0(4). (23) 

3) ^32 = 0(1). 

We briefly outline the calculation here: 

1) 2^i = 2%(2^ a + and 2^2% = (-2? + ^2^)2%. So 2^22?, - -S?i2 &? 2 = - 

22f2Sf a + ^(2Sf 2 2S^ - 25^2%) = Sf 2 Sf? - 22?2Sf 2 (By Lemma HTTl). 

Further computation also shows that, 

@21 = %2&* ~ 22f 2&? 2 = E/A.-EE'A./A* f °. ( 24 ) 

* * j 

Hence 7 ^ C£\T£ 2 , i.e., Jzff and =24 in general do not commute. Therefore, ^21 = 0 ( 1 ) 


( 22 ) 

= we will 


2) By direct computation we obtain that, 

C/d C/d C/d C/d C/d C/d I C/d C/d C/d _ C/dX C/)X c/d I c/d C/dX c/dX o r^.T O^.X 

^1 2 i - 242-^-1 - 2 >i — ^ 2-^1 

+ -4(2^2^42^2 + ^fT£fT£ 2 + SC 2 £Cf£Cf + - 22^42^4 - 22^42^4) (25) 

+ ^(25^2^2% + 2^'2Sf/’ - 225^2^2^) 

From Lemma Um we deduce that J£fJ£f££ 2 + ££ 2 J£f££f — 2£df£d 2 T£f = 0. The 0(^4) term 
can be simplified to < 2> 2 \T£f — T£fS> 2 \, but it is nonzero in general. 
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3) Following a similar calculation, we get, 


0 3 2 = &X&2&2 - 2.-^^ + J^^l 

= JSf^JgfaJgfa - 22S? 2 2S^2S? 2 + + \{£^X 2 X 2 - 2& 2 &?& 2 + J^J^) 

•~ 2 (26) 

_ c/) c/) OC^ C/) \ C/> C£ C/>X 

— JZsy ^2°^2 — £~Z'2~Z'\ ~&2 I 'Us2~£2~£\ 

= ^21^2-^2^21 = 0(1) 

Based on these explicit estimates, we now have for the original impulse method (J6J): 

Corollary II.3. When c\ — c 2 = ~,di = l,d 2 = 0, we have 

~ At 3 

&(ci, c 2 ; d u d 2 )~S = O(-) + ■■■ (27) 

£ z 


The important observation in this analysis is that the term 031 contains a large factor 
(dj). This motivates a different choice of the parameter: 


Corollary II.4. In the case when 



we have, 

&(ci, c-i i di, d?) — = 0(At 3 ) + • • • . 


(28) 


(29) 


In this case, we have abandoned the symmetry of the method, and chosen the parameters 
to eliminate the terms of the order G(At 2 ) and 0(^-) altogether. For better reference, we 
will call the original impulse method (J6]) the impulse I and the non-symmetric method (l28l) 

impulse II. 

In the case when k — 2, it is clear that the original impulse method is the only symmetric 
method. In order to obtain other symmetric methods, we need to consider splitting methods 
with more fractional steps. This will be discussed further in the next two sections. 


10 



2. Expansions of the error for k = 3 


When k — 3, we can choose d\ — (I 2 — d 3 — 0, C3 = ci, and C2 = 1 — 2c2 to form a 
symmetric integrator. In this case, the error 3% is given by, 


= f n/A-Xt 'A - ■va.xxa + a/a.A )^ 1 

+ (j - j + + i? 2 i? 2 ^)Ai 3 + 

= ( f -^ A(3 + (!-| + T 2 )s, - A(3 + "' 


(30) 


Based on the estimate in theorem III. 21 we choose to eliminate the term £^31, yielding, 


Corollary II.5. When C\ = c 3 = c 2 = d\ = d 2 = d 3 = 0, we have & — J? = 
^ 32 Af 3 + • • • . 


This method will be referred to as impulse III. For k = 3, there are also non-symmetric 
methods. But we will continue to consider the case k = 4. 


3. Expansions of the error for k = 4 

Finally, we will further explore the operator-splitting methods for the case k = 4. In this 
case, a symmetric method can be constructed by choosing ci, c 2 = | — Pl, c 3 = | — Ci, c 4 = 
ci; di, d 2 = 1 — 2o?i, d 3 = di, d± = 0. Up to At 4 terms, we have, 

^(ci, c 2 , c 3 , c 4 , di, d 2 , d 3 , d 4 ) - 

—(cidi — —di — cid 4 + — d 4 + — 2AfiAf2=Sfi + J^AfiAf^Af 3 

+ (^di — Cidi + — — )(Afi-Sf2«^2 — + Af^A^-SuJAt 3 

+ (^Cidi - - -cid? + -dl + —) 

x - Af 2 AfiAfiAfi + - 32^i2^ 2 A?i)At 4 (31) 

+ Mi - -d 3 - -dd? - -c?d! + -d? + —) 
v 1 1 g 2 1 i 2 i 1 4 1 lg ; 

x - Af) Af 2 + 2Af 1 ^ 2 ^i^2 - 2J^i^ 2 A?i)At 4 

A , 1 , 1 2 , 1 , 

v 2 8 2 1 48 ; 

x (j^j2? 2 2*? 2 j5f 2 - Af 2 ^ 2 ^ 2 ^i + 3^2^^^ - 32^i^ 2 A? 2 )Af 4 + • • • 
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Let us define, 


(3 A C/> C/> C/> C/> C/> C/> C/> C/> C/> C/> C/> O <V? c/> c/> c/> 

I ZC I ZC I ,X2 — ZC2~£\ZC\ZC\ t | ,X2^ l „Z. | — ■),X i ] | , 

042 =^ 2 ^\^\ - + 2jSf 1 j£f 2 j£f 1 jSf 2 - 2jgf 2 JgfiJgf 2 jgfi, (32) 

<G» A C/> C/> C/> C/> C/> C/> C/> C/> \ O C/> C/> C/> C/> n c/> c/> c/> c/> 

-^43 —,Z|.A'A-A ~~ ,Z 2 , : Z 2 ,Z2,Z| -f- •), : Z 2 „Z 2 ,Z|. : Z2 — •),Z 2 ,Zi,Z 2 . : Z 2 . 

A similar analysis yields, 

Theorem II.6. 77ie operators in (J32]) are o/ f/ie orders, 

1) ®41 = O(L), 

2) ®42 = 0(4), ( 33 ) 

£ z 

3) 043 = 0(1). 


In order to minimize the error, especially the first few terms of we will choose the 
parameters so that coefficients of 041 and 0 42 are zero (Notice that the coefficient of 031 is 
twice of the one of 0 41 , so it will be automatically zero). This leaves us with the nonlinear 
equations for Ci and dp. 


1 , 1 , 1 l9 1 l9 1 

~ c idi - -di - -ci^ + -dl + — - 0, 

, 3 , 1 , 2 1 2 , 1 , 2 1 

c i“i — g«i — 2 Cl “i — 2 < ' 1 ^ 1 + 4® 1 + Yq ~ 


(34) 


These equations can be simplified. In particular, d± is the only root of 6z 3 — 12z 2 + 6 z — 1, 
and Ci = ^di. 


Corollary II. 7. When ci = + -jy + |, d\ = 2c 1 , the coefficients of Q>, 31 , 03 2 , 04 i, 04 2 , 

and ^43 in CD are all zeros. 


The elimination of the 0 32 and 0 43 terms seem to be a coincidence. This method will 


1 


be referred to as impulse IV. With a direct calculation, one can show that Ci = -pr, 

2 ( 2 - 23 ) 

which surprisingly, coincides with the coefficients of the well known 4th order symplectic 


integrator 


15]. Of course, the symplectic method in 15] is based on the splitting of the 


kinetic and potential energy, while our splitting is between the fast and slow forces. 
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D. A Numerical Test: A nonlinear coupled oscillator 


Before we look further into the properties of the splitting methods, we present some 


numerical results for a nonlinear oscillator problem 


17], governed by the equations, 


-q + P\\6 - qf(9 -q) + -r (9 - q), 


0 = -p\\9-qf{9-q)--(9-q). 


Here q, 9 G M 2 . This is a Hamiltonian system with potential energy given by, 


(35) 


V(q, 6) = Ap - ,||* + |||» - ,11“ + i(||,|| - 1 ) 2 . 


(36) 


In our tests, the parameters are chosen as follows: e = 0.1, j3 = 0.1, q( 0) = (1,0), 9(0) = 
(1.01, 0), g(0) = (0,1), 0(0) = (0, 0.05), St = 0.01, and At = 0.12. 

First, we show the total energy computed from each method in Fig. [0 We find that 
the new impulse methods have much better performance in the energy conservation: The 
fluctuation is much smaller than the original impulse method (impulse I). For problems 
where the energy is more relevant than the actually trajectories, e.g., producing various 
statistical ensembles, the new methods seem to be more promising. Among the new impulse 
methods, the method III seems to have the best results. It is clear, however, much deeper 
analysis is needed to understand the accuracy of the methods toward computing different 
quantities. We will present some preliminary analysis in the next section. 

In Fig. H we show the error of q± for the numerical approximations obtained from the 
impulse methods. The error is estimated by comparing the approximate solutions to a 
solution computed with very small step size. We observe that the accuracy is gradually 
improved for the impulse methods I to IV. However, the error would grow in all cases, 
which can be attributed to the Lyapunov instability inherent in most Hamiltonian systems. 

In Fig. [3l we show the momentum computed from the four methods. Surprisingly, they 
exhibit similar accuracy, and the new impulse methods show little improvement. 
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FIG. 1. Comparison of the energy conservation for the four methods. 
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FIG. 2. Comparison of the error in q\ for the four methods. 
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FIG. 3. Comparison of the error in p± for the four methods. 


E. Some preliminary analysis 


1. Energy conservation of multiscale Hamiltonian systems 


For those ODEs (J2j) that come from Hamiltonian systems, i.e., / = —M X VW and 
g = —M -1 V V, we define two Hamiltonians that correspond to the splitting of the operator 

H 2 =W{x), 


H i =— V(x) + -p 1 AT V 


(37) 


e* ' ■ T 

Here M is the mass matrix, and p = Mv is the momentum variable. 

Due to the approximation, the energy associated with the dynamical system will not be 
exactly conserved. However, one of the celebrated results in geometric integrators is that an 
approximate Hamiltonian often exists, and it is conserved exactly by the numerical method 


( 121 , 


17-20], 


For impulse I the energy conservation property can be analyzed using the backward 
analysis |12l . 20], which asserts that the approximate solution is a more accurate solution 
of another Hamiltonian system with a Hamiltonian H $, known as the shadow Hamiltonian. 
For the impulse method I, the analysis shows that the shadow Hamiltonian, up to the order 
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At 2 , is given by [12], 


At 2 At 2 

H s = H 1 + H 2 + — {{H 2 , H,}, H,} - — {# 2 , {H 2 , H,}}. 


( 38 ) 


This is because the operator approximation can be written as, 


e iA*jS? 2e AtS?i e §Atj2?2 _ e 2’At+^ r 2>3i-^t2>32+~ 


(39) 


Here, = [[Jf 2 , Afi] and ^ 32 = [Af 2 , [-& 2 , ]]; [ ] and { } stand for the commutator 

(Lie derivative) and Poisson bracket, respectively. 

In particular, we have that, 


{H 2 ,H 1 } = -p 1 M~ l VW(x), 

{{H u H 2 },H 2 } =p T M- 1 V 2 WM~ 1 p- ^VW t M~ 1 VV(x), 
=VW t M~ 1 VW. 


(40) 


As a result, the conservation of the energy at this level is dominated by the term, 

which is due to the presence of fi? 31 in the error 

In contrast, the same calculation for the impulse method II yields, 

1 7 At 2 

H S = H+ H 2 , {H 2 , H,}} = H + 0(At 2 ) +■■■ (41) 

As a result, the better energy conservation can be attributed to the elimination of the ^31 
term in the error. 


2. Resonance instability 


Another outstanding issue raised by previous works is the resonance, which occurs for 
certain choices of the slow time step At { 3 , 6|. Following the analysis in [6, 8|, and in 
particular the example in Q, we consider a scalar problem where f(x) = — (-|) 2 x and 
g(x) = —7t 2 x. In Fig. HI we show the spectral radius of the propagation matrix. We observe 
that the new impulse methods exhibit similar resonance phenomena: When the large step 
size At is around an integer multiple of half of the period (T=2) associated with the fast 
scale, instability occurs. 
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FIG. 4. Spectral radius of the propagation matrix. 


III. ANOTHER EXAMPLE: DYNAMICS OF OCTANE 


Here, we consider the Octane molecule with 26 atoms. The impulse methods have been 
implemented within TINKER { 4 ]. To properly quantify the error, the ‘exact’ solution is 
represented by the solution computed with the Verlet’s method with small step size 10 ~ 5 ps. 
In the impulse method, we choose At = 2.4 x 10 -4 ps and 5t = At/24. All the simulations 
are conducted for 6 ps period. 

In Fig. 0 we show the total energy computed from impulse methods I to III. Again we 
observe that the new impulse methods have much less fluctuation of the energy, indicating 
a better energy conservation property. 

Next we look at the error in the position of the first atom (first component). The results 
are shown in Fig. [6] There are some improvement of the accuracy from impulse methods I 
to impulse III. Such improvement is also observed in the velocity (first component iq), as 
can be seen in Fig. [3 


IV. SUMMARY AND DISCUSSION 

We have developed some new impulse methods for the numerical approximation of molec¬ 
ular systems with multiple time scales, which are represented by interactions of different 
magnitude. Motivated by the operator-splitting approach of the original impulse method 


17 






























Impulse I: Total energy 
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FIG. 5. Comparison of the energy conservation. 
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FIG. 6. Comparison of the error in x\. 


, 2], we sought general splitting methods that involve more fractional steps. A novel aspect 
in our approach is the systematic procedure for finding an expansion of the error, which in 
turn sheds light on the selection of the coefficients so that the accuracy can be improved. 

For multiscale ODEs, our analysis revealed that the terms in the error can depend on 
both the large time step At, and s, which represents the separation of the scales. Based on 
the order of the first few terms, we choose the coefficients so that the terms with the largest 
magnitude are eliminated. This leads to several splitting methods that can be viewed as 
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Impulse I: Error in velocity 
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FIG. 7. Comparison of the error in v\. 


generalized impulse methods. Numerical tests have been conducted, which have confirmed 
the improved accuracy. The biggest improvement has been observed in the energy conserva¬ 
tion. This has been analyzed with some preliminary study of the modified Hamiltonian. It is 
thus expected that the new methods would produce better results when the micro-canonical 
ensemble distribution is of interest. 


This approach, however, is by no means complete. From a practical viewpoint, the 
Hamiltonian system (J2]) models a system in isolation. In practice, often of interest are 
extended systems, where the external conditions are modeled by introducing additional 
variables, e.g., heat and pressure bath, or by introducing stochastic forces, e.g., the Langevin 
dynamics. Operator-splitting methods have been widely used for extended systems, e.g., in 
1.4], and the impulse methods have also been applied to Langevin dynamics 21] as well. 
Extending the current framework to those problems might produce new integrators with 
other capabilities, and it will be explored in our future works. 

The problem considered in this work belongs to stiff ODEs, for which many numeri¬ 
cal methods have been developed, e.g., implicit Runge-Kutta methods and BDF methods, 
and they can be found in standard textbooks [if]]. Meanwhile, there have been significant 


recent progress in developing efficient computational methods for such dynamica 
with multiple time scales, e.g., the heterogeneous multiscale method (HMM) 


equation-free method 


systems 

425]. the 

26], the FLAVOR method 27], the reversible averaging integrator 
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[171 ] . etc. These methods demonstrate resemblance to the impulse method in that they in¬ 
troduce multiple time steps (St and At) to capture the multiple scales. On the other hand, 
an averaging procedure is usually involved on the fastest time scale to compute an effective 
force on quantities that evolve on the slow time scale. In addition, some of these methods 
assume the existence and explicit form of slow variables. At this point, we are not aware of 
the application of these methods to biomolecular modes. 
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